\defmodule {CIRProcess}

This class represents a \emph{CIR} (Cox, Ingersoll, Ross) process \cite{fCOX85a}
 $\{X(t) : t \geq 0 \}$, sampled at times $0 = t_0 < t_1 < \cdots < t_d$.
This process obeys the stochastic differential equation
\begin{equation}
   dX(t) = \alpha(b - X(t)) dt + \sigma\sqrt{X(t)}\, dB(t)
                                               \label{eq:cir}
\end{equation}
with initial condition $X(0)= x_0$,
where $\alpha$, $b$ and $\sigma$ are positive constants,
and $\{B(t),\, t\ge 0\}$ is a standard Brownian motion
(with drift 0 and volatility 1).
This process is \emph{mean-reverting} in the sense that it always tends to
drift toward its general mean $b$.
The process is generated using the sequential
technique \cite[p. 122]{fGLA04a}
\begin{latexonly}
\begin{equation}
   X(t_j) = \frac{\sigma^2\left(1 - e^{-\alpha(t_j - t_{j-1})}\right)}{4\alpha}
   \chi^{\prime\,2}_\nu\left(\frac{4\alpha e^{-\alpha(t_j - t_{j-1}) } X(t_{j-1})}
   {\sigma^2\left(1 - e^{-\alpha(t_j - t_{j-1})}\right)}\right),
                                    \label{eq:cir-seq}
\end{equation}
\end{latexonly}
\begin{htmlonly}
\begin{equation}
   X(t_j) = \left\{}{\sigma^2\left(1 - e^{-\alpha(t_j - t_{j-1})}\right)/ {4\alpha}\right\}
   \chi^{\prime2}_\nu\left(\{4\alpha e^{-\alpha(t_j - t_{j-1}) } X(t_{j-1})\}/
   \{}{\sigma^2\left(1 - e^{-\alpha(t_j - t_{j-1})}\right)}\}\right),
                                    \label{eq:cir-seq}
\end{equation}
\end{htmlonly}
where $\nu = 4b\alpha/\sigma^2$, and $\chi^{\prime\,2}_\nu(\lambda)$ is a noncentral
chi-square random variable with $\nu$ degrees of freedom and noncentrality
parameter $\lambda$.

\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        CIRProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;

\end{hide}

public class CIRProcess extends StochasticProcess \begin{hide} {
    private RandomStream stream;
    protected ChiSquareNoncentralGen gen;
    protected double alpha,
                     beta,
                     sigma,
                     nu;     // number of degrees of freedom
    // Precomputed values
    protected double[] parc,
                       parlam;
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public CIRProcess (double x0, double alpha, double b, double sigma,
                      RandomStream stream) \begin{hide} {
      this (x0, alpha, b, sigma, new ChiSquareNoncentralGen (stream, null));
    }\end{hide}
\end{code}
\begin{tabb} Constructs a new \texttt{CIRProcess} with parameters
$\alpha = \texttt{alpha}$, $b$, $\sigma = \texttt{sigma}$ and initial value
$X(t_{0}) = \texttt{x0}$. The noncentral chi-square variates
  $\chi^{\prime2}_\nu(\lambda)$ will be
generated by inversion using the stream \texttt{stream}.
\end{tabb}
\begin{code}

   public CIRProcess (double x0, double alpha, double b, double sigma,
                      ChiSquareNoncentralGen gen) \begin{hide} {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      nu = 4.0*b*alpha/(sigma*sigma);
      this.gen   = gen;
      stream = gen.getStream();
    }\end{hide}
\end{code}
\begin{tabb} The noncentral chi-square variate generator  \texttt{gen}
 is specified directly instead of specifying the stream.
 \texttt{gen} can use a method other than inversion.
\end{tabb}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}\begin{hide}

   public double nextObservation() {
      double xOld = path[observationIndex];
      double lambda = xOld * parlam[observationIndex];
      double x;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = parc[observationIndex] *
             ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      observationIndex++;
      path[observationIndex] = x;
      return x;
   }

\end{hide}
   public double nextObservation (double nextTime) \begin{hide} {
      double previousTime = t[observationIndex];
      double xOld = path[observationIndex];
      observationIndex++;
      t[observationIndex] = nextTime;
      double dt = nextTime - previousTime;
      double x = nextObservation (xOld, dt);
      path[observationIndex] = x;
      return x;
    }\end{hide}
\end{code}
\begin{tabb} Generates and returns the next observation at time $t_{j+1} =
 \texttt{nextTime}$, using the previous observation time $t_{j}$ defined earlier
(either by this method or by \texttt{setObservation\-Times}),
as well as the value of the previous observation $X(t_j)$.
\emph{Warning}: This method will reset the observations time $t_{j+1}$
for this process to \texttt{nextTime}. The user must make sure that
the $t_{j+1}$ supplied is $\geq t_{j}$.
\end{tabb}
\begin{code}

   public double nextObservation (double x, double dt) \begin{hide} {
      double c = -Math.expm1(-alpha * dt) * sigma * sigma / (4.0*alpha);
      double lambda = x * Math.exp(-alpha * dt) / c;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = c * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = c * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = c * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      return x;
    }\end{hide}
\end{code}
\begin{tabb} Generates an observation of the process in \texttt{dt} time units,
assuming that the process has value $x$ at the current time.
Uses the process parameters specified in the constructor.
Note that this method does not affect the sample path of the process
stored internally (if any).
\end{tabb}
\begin{hide}\begin{code}

   public double[] generatePath() {
      double xOld = x0;
      double x, lambda;
      int j;

      if (gen.getClass() == ChiSquareNoncentralPoisGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else if (gen.getClass() == ChiSquareNoncentralGamGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }
      }

      observationIndex = d;
      return path;
   }

   public double[] generatePath (RandomStream stream) {
      gen.setStream (stream);
      return generatePath();
   }
\end{code}
\begin{tabb} Generates a sample path of the process at all observation times,
 which are provided in array \texttt{t}.
 Note that \texttt{t[0]} should be the observation time of \texttt{x0},
 the initial value of the process, and \texttt{t[]} should have at least $d+1$
 elements (see the \texttt{setObservationTimes} method).
\end{tabb}\end{hide}
\begin{code}

   public void setParams (double x0, double alpha, double b, double sigma) \begin{hide} {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      nu = 4.0*b*alpha/(sigma*sigma);
      if (observationTimesSet) init(); // Otherwise not needed.
    }\end{hide}
\end{code}
\begin{tabb}
Resets the parameters $X(t_{0}) = \texttt{x0}$, $\alpha = \texttt{alpha}$,
 $b = \texttt{b}$ and $\sigma = \texttt{sigma}$ of the process.
\emph{Warning}: This method will recompute some quantities stored internally,
which may be slow if called too frequently.
\end{tabb}
\begin{code}

   public void setStream (RandomStream stream) \begin{hide} { gen.setStream (stream); }\end{hide}
\end{code}
\begin{tabb}
Resets the random stream of the noncentral chi-square generator to \texttt{stream}.
\end{tabb}
\begin{code}

   public RandomStream getStream() \begin{hide} { return gen.getStream (); }\end{hide}
\end{code}
\begin{tabb}
Returns the random stream of the noncentral chi-square generator.
\end{tabb}
\begin{code}

   public double getAlpha() \begin{hide} { return alpha; }\end{hide}
\end{code}
\begin{tabb}
Returns the value of $\alpha$.
\end{tabb}
\begin{code}

   public double getB() \begin{hide} { return beta; }\end{hide}
\end{code}
\begin{tabb}
Returns the value of $b$.
\end{tabb}
\begin{code}

   public double getSigma() \begin{hide} { return sigma; }\end{hide}
\end{code}
\begin{tabb}
Returns the value of $\sigma$.
\end{tabb}
\begin{code}

   public ChiSquareNoncentralGen getGen() \begin{hide} { return gen; }\end{hide}
\end{code}
\begin{tabb}
Returns the noncentral chi-square random variate generator used.
The \texttt{RandomStream} used for that generator can be changed via
\texttt{getGen().setStream(stream)}, for example.
\end{tabb}
\begin{code} \begin{hide}

   protected void initArrays(int d) {
      double dt, c;
      for (int j = 0; j < d; j++) {
         dt = t[j+1] - t[j];
         c = -Math.expm1(-alpha * dt)*sigma*sigma/(4.0*alpha);
         parc[j] = c;
         parlam[j] = Math.exp(-alpha * dt) / c;
      }
   }

   // This is called by setObservationTimes to precompute constants
   // in order to speed up the path generation.
   protected void init() {
       super.init();
       parc = new double[d];
       parlam = new double[d];
       initArrays(d);
    }\end{hide}
\end{code}
\begin{code}\begin{hide}
} \end{hide}
\end{code}
